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Abstract 

Extracorporeal Shock Wave Therapy (ESWT) is a noninvasive treatment for a variety of 
musculoskeletal ailments. A shock wave is generated in water and then focused using an acous- 
tic lens or reflector so the energy of the wave is concentrated in a small treatment region where 
mechanical stimulation enhances healing. In this work we have computationally investigated 
shock wave propagation in ESWT by solving a Lagrangian form of the isentropic Euler equa- 
tions in the fluid and linear elasticity in the bone using high-resolution finite volume methods. 
We solve a full three-dimensional system of equations and use adaptive mesh refinement to 
concentrate grid cells near the propagating shock. We can model complex bone geometries, 
the reflection and mode conversion at interfaces, and the the propagation of the resulting shear 
stresses generated within the bone. We discuss the validity of our simplified model and present 
results validating this approach. 

1 Introduction 

Extracorporeal shock wave therapy (ESWT) is a noninvasive treatment for musculoskeletal condi- 
tions such as bone fractures that fail to heal (non-unions), necrotic wounds, and strained tendons [61, 
45] . In this treatment a shock wave is generated in water and then focused using an acoustic lens or 
reflector so that the energy of the wave is concentrated in a small treatment region. This technique 
has been used since the 1980's, more widely in Europe and Asia than in the US, where it is still 
considered experimental and has limited FDA approval. 

Although the underlying biological mechanisms are not well understood [47], the mechanical 
compressional and/or shear stress caused by the propagating shock wave is thought to stimulate 
healing [47, 60, 29, 44, 58, 48, 49, 18, 13, 30, 11, 26]. There have been a variety of experimental 
studies to investigate the effects of extracorporeal shock waves on urinary stones of varying hardness, 
as well as other hard tissues such as bone and its surrounding tissues (cartilage, tendon and fascia). 
In particular, it has been shown that high extracorporeal shock wave energy actually fractures rat 
bones, but lower applied energy levels stimulated osteogenesis [32, 59, 63] We should also note that 
in a recent review article, Zelle, et. al. concluded that the majority of clinical studies were done at 
too high a level to enable a clinical recommendation [64]. However, they did conclude that shock 
wave therapy seems to stimulate the healing process in delayed unions and nonunions [3]. There 
are a number of other biological mechanisms in addition to stress that potentially play a role in the 
body's response to ESWT. The focus of this study, however, is on mechanical stress deposition and 
computational tools for studying this aspect. 

The medical shock wave devices are similar to those used for extracorporeal shock wave lithotripsy 
(ESWL), a widely- used non-surgical treatment for kidney stones in which the focused shock waves 
have sufficient ampHtude to pulverize the kidney stone. In shock wave therapy the amplitudes 
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a) b) 

Figure 1: Cartoon of the Dornier HM3 Lithotripter. In a) the spherical wave is generated at Fl, 
reflects off the ellipsoid and the reflected wave focuses at F2. In b) the diagram illustrates the 
creation of the edge waves at the corner of the ellipsoid and the contribution of negative pressure to 
the tail of the ESWT pressure wave. 

are generally smaller and the goal is mechanical stimulation rather than destruction, although in 
some applications such as the treatment of heterotopic ossifications (HO) (see Section 5.4) larger 
amplitudes may be used. 

Figure 1 shows the geometry of a laboratory shock wave device modeled on the clinical Dornier 
HM3 lithotripter. The three-dimensional axisymmetric geometry consists of an ellipsoidal reflector 
made out of metal and a cavity filled with water. A spark plug at the focus of the ellipse marked 
Fl generates a bubble which collapses and creates a spherical shock wave that reflects and focuses 
at F2. The major and minor axes of the ellipsoid in the HM3 are a = 140mm and b = 79.8mm, 
respectively. The foci of this ellipse are at (±115, 0, 0) and the reflector is truncated at 100mm from 
Fl, or (-10,0,0). 

In the laboratory, this reflector is immersed in a bath of water and objects can be placed at the 
second focus of the ellipsoid, F2. This device is in use at the Center for Industrial and Medical 
Ultrasound (CIMU) at the University of Washington Applied Physics Laboratory and we have used 
this geometry in order to compare directly with some laboratory experiments. Some preliminary 
comparisons were presented in [22]. 

Computationally, we use this geometry to calculate the initial condition by solving two-dimensional 
axisymmetric Euler equations with the Tammann equation of state (see Section 2). These initial 
conditions are then fed into a full three-dimensional calculation near the focus at F2. 

In addition to the RMS, we have also used the geometry of the hand- held Sanywave device used 
in clinical studies by our collaborator Dr. Michael Chang. Some sample calculations related to the 
study of HOs are presented in Section 5.4. 

In each case, the ESWT pressure wave form that is generated has a similar shape. There is a 
sharp increase in pressure from atmospheric pressure (~ O.lMPa) to a peak pressure ranging from 
35 to 100 MPa over a very short rise time (~ 10 ns), followed by a decrease in pressure to ~ — 10 
MPa over ~ 5/is. The negative fluid pressure in the tail can lead to cavitation bubbles, as discussed 
below. 

Computational models for shock wave propagation and focusing can aid in the study of ESWT. 
In particular, there are many open questions concerning the interaction of shock waves with complex 
three-dimensional geometries such as bone embedded in tissue. 

Because of the difference in material properties, a wave hitting the tissue/bone interface will 
be partially reflected, and the transmitted wave will have a modified strength and direction of 
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propagation. This can greatly affect the location and size of the focal region as well as the peak 
pressure amplitude. Moreover, although the shock wave is primarily a pressure wave in soft tissue 
(which has a very small shear modulus), at a bone interface mode conversion takes place and shear 
waves as well as compressional waves are transmitted into the bone, generating a dynamically applied 
load. 

Bone healing is thought to be regulated in part by mechanical factors [44, 29, 58, 50, 29, 26, 58]. 
Several studies have shown that the application of cyclic compressive and shear displacements can 
enhance healing through increased callus formation and ossification [44, 50, 51, 55, 48, 62]. The 
results also indicate that treatment is also dependent upon the rate, mode and magnitude of the 
stress deposition [44], as well as the gap size [14]. 

Carter, et. al. [11], as well as, and Claes and Heigele [13], proposed a model for skeletal 
tissue development based on hydrostatic pressure and tensile displacements [13]. Other research has 
proposed a different model for skeletal tissue formation based on shear strain and fluid flow [49, 18]. 
Augat, et. al. [1], found that tensile displacements are not effective in enhancing bone formation. 
This was further validated whe Isaksson, et. al. [30] investigated the models in [13, 11, 49, 18] and 
found that shear strain and fluid flow, were more accurate predictors of bone growth. However, no 
single model was able to predict certain features of the bone formation and healing process [44], 
highlighting the need for further research in this area. 

The shear waves generated at the fluid/solid interface have also been shown to be important in 
the effective break up of kidney stones [54, 23]. An additional effect of ESWT is the formation and 
collapse of cavitation bubbles that can cause tissue damage. While the shock wave is a compression 
wave, it is followed by a rarefaction wave of expansion, and in the tail the fluid pressure typically 
drops to negative values. Reflection at interfaces can lead to enhanced regions of expansion and to 
sufficiently negative pressures that cavitation bubbles can form [43, 56, 17]. 

To better understand all of these effects, it is desirable to have a three-dimensional computational 
model that can simulate the focusing of nonlinear shock waves and their interaction with arbitrarily 
complex interfaces between different materials. 

In this paper we present an approach to this problem that has allowed the study of some of these 
issues in a simplified context. In particular, we consider an idealized situation in which soft tissue is 
replaced by water, ignoring its viscoelastic properties, and modeled by the nonlinear compressible 
Euler equations with the Tammann or Tait equation of state. This has been used for prior ESWT 
work in water as well as biological- like materials [31, 46]. Bone is modeled as an isotropic and 
homogeneous linear elastic material [24, 33]. 

In reality, soft tissue and bone are very complex multiscale materials with microstructures, inho- 
mogeneities, and anisotropic properties. Any attempt to model the biological effect of shock wave 
propagation through such materials may require a more sophisticated and detailed model than used 
here. However, we believe that many of the macro-scale shock propagation issues discussed above 
can be adequately and most efficiently studied with a simplified model of the form considered here, 
since the dominant effect we hope to capture is the reflection and transmission of waves at interfaces 
between materials. 

The compressible Euler equations with the Tammann equation of state (see Section 2.1) in two- 
dimensional axisymmetric geometry is used to model the initial formation of the focusing shock 
wave. These initial conditions are then fed into a code that uses a simpler nonlinear model, the Tait 
equation of state, in a three-dimensional simulation of the fluid. The compressible fluid equations 
are written using a Lagrangian formulation that easily couples to the isotropic linear elasticity 
equations used in the bone-like material. The resulting equations have the same form everywhere, 
with a different stress-strain relationship in the different materials. 

A high- resolution finite volume method is used to solve these equations. We use the wave- 
propagation algorithms described in [40] and implemented in Clawpack [39]. These are Godunov- 
type methods for the hyperbolic system that use solutions to the Riemann problem between adjacent 
grid cells to determine a set of waves used to update the solution, and second-order correction terms 
with slope limiters are added to resolve the nearly discontinuous shock waves with minimal smearing 
or nonphysical oscillation. 
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These methods are used on a purely rectangular Cartesian grid. Each grid cell has associated 
with it a set of material parameters determining the material in the cell, in a unified manner so 
that both fluid and solid can be modeled. Complex geometry is handled by using appropriate 
averaged values of these parameters in cells that are cut by the interface. This is described further 
in Section 4.4. Averaging across the interface works quite well when the material properties are 
sufliciently similar and in Section 4.4 we show that this is the case even for fluid/solid boundaries 
of the type we consider. 

We also use patch-based adaptive mesh refinement (AMR) to concentrate grid cells in regions 
where they are most needed to resolve features of interest. The Clawpack software contains AMR 
software in both two and three space dimensions and this software has been used directly for the 
two-dimensional axisymmetric computations of the initial shock wave described in Section 5. For 
the three-dimensional problem we have used ChomboClaw [10], an interface between Clawpack and 
the Chombo code developed at LBL [35], which provides an implementation of AMR on parallel 
machines using MPI. Using ChomboClaw, the code originally developed using Clawpack was easily 
converted into a code that was run on an NSF TeraGrid machine at Texas Advanced Computing 
Center (TACC) and tested using up to 128 processors. 

Extensive laboratory experiments have been performed on shock wave devices to measure the 
wave form of shock waves produced by various devices, the shape of the focal region, the peak ampli- 
tudes of pressure observed in these regions, and other related quantities. Most of these experiments 
have been done in a water tank where the shock wave propagates and focuses in a homogeneous 
medium where measurements are easily done, or with phantoms (acrylic objects with well under- 
stood photoelastic properties) that are placed in the water as a proxy for bones or kidney stones, 
with instrumentation such as pressure gauges or photographs used to explore the interaction of the 
shock wave with the object. In some cases high-speed photographs of the shock wave have been 
obtained. Creating phantoms from clear bi-refringent materials and using polarized light it is even 
possible to photograph the shock wave propagating through the object [53]. We have used some of 
these experiments to help validate our numerical approach [22]. 

Other researchers have also developed computational models for shock wave therapy and lithotripsy. 
In prior work the pressure field has been modeled using linear and nonlinear acoustics as well as the 
Euler equations with the Tait equation of state. Hamilton [27] used linear geometrical acoustics, 
which holds under the assumption of weak shock strength, to calculate the reflection of the spherical 
wave. The diffraction of the wave at the corner of the reflector was calculated using the Kirchoff 
integral method. Christopher's [12] model of the HM3 lithotripter used Hamilton's result as a start- 
ing point and considered non planar sources. Coleman et al [17], Averkiou and Cleveland [2] used 
models based on the KZK equation. Tanguay [56] solved the full Euler equations and incorporated 
cavitation effects as well as the edge wave. 

Our approach differs from these in that we consider the wave propagation in both the fluid and 
solid by solving a single set of equations that can model both materials. This approach allows us 
to investigate not only compression and tension effects of ESWT, but also the propagation of shear 
waves in the solid. Sapozhnikov and Cleveland [15] have investigated the effect of shear waves on 
spherical and cylindrical stones using linear elasticity with a plane wave initial condition. This initial 
condition is an unfocused wave, which yields good results for small objects, but would fail to capture 
the full ESWT pressure wave interaction with three-dimensional bone geometries. 

2 Model equations 

To accurately model shock wave formation and propagation it is generally necessary to use nonlinear 
equations of compressible flow. In this work we use nonlinear equations for compressible liquids in 
the fluid domain (water or soft tissue) and linear elasticity in the solid domain (bone). The nonlinear 
compressible equations are written in a Lagrangian framework in terms of a reference conflguration, 
as is done for the linear elasticity equations. This allows both sets of equations to be written 
in the same form. We apply flnite volume methods to this form of the equations so that a single 
computational grid (or set of nested grids with AMR) can be used over the entire domain. Interfaces 
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between fluid and solid are represented by choosing averaged material parameters in each grid cell, 
as discussed further in Section 4.4. 

The system of equations we solve has the general form of a hyperbolic system of 9 equations 



Qt + /(g, X, y, z)^ + g{q, x, y, z)y + h{q, x, y, z)^ = 0, 



(1) 



where the vector q consists of the 6 components of the symmetric strain tensor followed by the 
momenta, and the fluxes in general may be spatially varying based on material properties: 



^11 



^33 
^23 

pu 
pv 
pw 



f{q,x,y,z) 



u 



v/2 


w/2 
^11 



(J 



g{q,x,y,z) 





V 



u/2 
w/2 


^23 



h{q,x,y,z) 






w 

v/2 

u/2 

^33 



(2) 



In these expressions, p = p{x^y^z) is the density of the material (the "background density" inde- 
pendent of the wave propagating through the material) and the stress tensor a = a{q, x, y, z) is in 
general a spatially varying function of q^ linear in the solid and nonlinear in the fluid. 

Within the fluid domain <j = — p/, where p is the scalar pressure and / is the identity matrix. 
The pressure is a nonlinear function of the strain as discussed further below. In the solid domain, 
cr is a linear function of e and is non-diagonal, allowing us to model the propagation of shear waves 
as well as compressional waves. 

In Section 2.1 below we present the compressible fluid equations in their standard Eulerian form 
(the Euler equations) and discuss two possible equations of state, the Tammann EOS and the simpler 
Tait EOS in which the pressure is a function of density (or strain) alone, allowing us to drop the 
energy equation from the Euler equations. Then in Section 2.2 we rewrite these equations in the 
Lagrangian form given above. This can be done when modeling ESWT because the deformations 
are sufliciently small that the geometric nonlinearity of the equations can be ignored, adopting a 
Lagrangian frame and only considering the nonlinearity of the stress-strain relation as given by the 
equation of state. 

In Section 2.3 we discuss the linear elasticity model use to model bone. 



2.1 Compressible fluids in Eulerian form 

Much of the previous work on ESWT has been centered around the use of the Euler equations with 
the Tait or Tammann equations of state. These equations of state are typically used for modeling 
underwater explosions like the spark plug source of the lithotripter device [27, 31]. In this section we 
discuss the full Euler equations and proceed to show why the Tait equation of state is sufficient for 
modeling ESWT. Since this equation of state is a function only of the density, and can be rewritten 
as a function of strain, we show in Section 2.2 how it can be modeled within the framework of 
elasticity, which enables us to model both the fluid and solid with the single system of equations 
given above. 

In three space dimensions the Euler equations take the form 
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Several of the problems we investigated are axially symmetric and this enabled us to reduce 
the three-dimensional equations to a two-dimensional form. If we first rewrite the equations in 
cylindrical coordinates (r, 0, z) and assume no variation and zero velocity in the ^-direction, the 
system we obtain is reduced to two variables, r and z. The equations are 



(4) 



where and Wz denote the velocities in the r and z directions, respectively. These equations are 
of the same form as the two-dimensional Euler equations, with the addition of geometric source 
terms that are a result of the variable transformation. The source terms are never evaluate at r = 
since we are using a finite volume method where quantities are evaluated at cell-centers, that is, the 
smallest value of r in a calculation is Ax/2. We prefer to keep the equations in conservation form, 
so they can be eflSciently solved using finite volume methods. 

In order to solve the system (3) or (4), we need to close the system with a relation between the 
pressure and conserved variables. The Tammann EOS [31] is applicable to a wide range of liquids, 
even with very strong shock waves. This equation of state has the form 



P = e) = (7 - l)pe - 



(5) 



where p, p and e are the pressure, density and specific internal energy, respectively, while 7 and Poo 
are constants depending on the fiuid. If Poo = this is the standard EOS for an ideal gas, with 7 
generally satisfying 1 < 7 < 5/3, while for water 7 « 7.15 and Poo ~ 300MPa. For sufficiently weak 
shocks, this can be approximated by the Tait equation of state. 



p = p{p) = B 
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(6) 



where 5 is a pressure term that is a weak function of entropy, but is typically treated as a constant, 
and corresponds to from (5) while n corresponds to 7. In our work we take B = 300MPa and 
n = 7.15. 

It has been common practice to use the Tait EOS in shock wave therapy and lithotripsy models 
[52, 46]. This has been justified by noting studies that show that entropy changes across the shock 
are very small even up to pressure jumps of 200 MPa [46], which is beyond the range used in ESWT. 
To verify this assumption, we performed computational experiments to compare the Tammann and 
Tait equations 

of state for typical ESWT shock waves. Since we have used the f-wave approach in our compu- 
tational model, we can solve 4 with a spatially-varying equation of state. We set up an experiment 
where the resulting shockwave (generated using the Tammann equation of state), was over 150 MPa. 
Figure 2 a) shows the results from this experiment. The black solid curve is the result from solving 
with the Tammann EOS in the entire domain. The blue dashed curve shows the result gotten by 
switching to the Tait EOS at x = 50. This enabled us to compare the two equations of state with the 
exact same initial condition. There is a small disagreement at x = 50 caused by a slight refiection at 
the interface due to the change in the equation of state. Otherwise, the pressure profiles are nearly 
identical, giving confidence that the calculations we are interested in can be done by solving the 
Euler equations with the Tait equation of state. This allows us to drop the equation for energy and 
obtain the simplified system 
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Zoomed in Pressure Comparison 




Lagrangian vs. Eulerian Calculation 




- Lagrangian 
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Figure 2: a) Comparison of a pressure wave calculation performed using both the Tait (blue dashed 
curve) and Tammann (black curve) equations of state. The results are nearly identical, b) Com- 
parison of the pressure pulse at F2 obtained in the Euler calculation (blue dashed curve) and the 
Lagrangian calculation (black curve). It is clear that the two sets of equations give good agreement. 
The wave in the Lagrangian case is slightly attenuated, but this may be due to error in initializing 
the calculation. In these calculations Ax = 0.5mm. 



2.2 Compressible fluids in Lagrangian form 

In the case of a fluid where the shear modulus is zero, the stress tensor can be written as (j{e) = — p/, 
where p is the pressure in the fluid and / is the identity matrix. In the case of ESWT, the pressure 
only depends on changes in the density, and we can write p{e) as a function of the strain tensor e. 
Consider the movement of a material with respect to a reference configuration and let 6 = ((5^, (5^, 5^) 
be the infinitesmal displacement. 

In three space dimensions, the full strain tensor is 

I) hisi + s^) \ 

(8) 




where subscripts denote partial derivatives. 

In the case of small deformations, we have that 
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where po is the equilibrium density. 

If we insert this into the Tait equation of state (6) we get 
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(10) 



Using the Lagrangian form is only valid in the case where the displacements are small, so we 
calculated the maximum value of the displacements in a two-dimensional axisymmetric calculation 
with the Euler equations. We found that for the maximum peak pressure of 50MPa, the correspond- 
ing maximum velocity was 10~^m/s. We then calculated the maximum displacement by integrating 
the velocity over the time of the calculation and found this to be on the order of 10~^mm. The size 
of the grid cell is on the order of 10~^mm, so the displacements are 4 orders of magnitude smaller 
than the width of the grid cells. It is therefore reasonable to assume that the density in each grid 
cell is essentially constant and that the Lagrangian framework of the elasticity equations will be 
valid for the fluid. 

To test this, we took the same initial condition for the 2D axisymmetric Euler equations with 
the Tammann equation of state and the corresponding 2D axisymmetric Lagrangian form of the 
equations with the Tait equation of state and measured the pressure at the focus, F2. The results in 
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Figure 3: This series of figures shows a pressure gauge measurement at F2 of different versions of 
the Tait EOS at different amphtudes. The triangular markers indicate the fuh nonhnear Tait EOS, 
the sohd hne is a hnearized version and the square markers are a quadratic version. The hnearized 
versions of the EOS work reasonably well at small amplitudes, but it's clear from figure c) that as 
the pressures increase to those observed in ESWT, the full nonlinear equation of state must be used. 



Figure 2 b) demonstrate reasonably good agreement between the two cases, but the Lagrangian form 
is slightly attenuated. This may be due to conversion of the initial condition from the conserved 
variables in the Euler equations (4) to those in the elasticity equations (1). 

Since the displacements are small, we also considered the possibility that nonlinearity in the fluid 
could be ignored, so we could instead using a linearized version of the Tait equation of state. Then 
we would be able to simply use the linear elasticity equations throughout the domain, in both the 
fluid and solid materials. If we assume a small perturbation to the strain, e + ^e, we can expand the 
Tait EOS (6) as a Taylor series about e, 

p{e + 5e)=po+p'{5e)e+^^e^ + --- (11) 

If we keep the first two terms of the expansion, the EOS has been linearized and we will call this 
the "linear Tait EOS." Similarly, we will refer to the equation obtained by keeping the first three 
terms of the expansion as the "quadratic Tait EOS." One-dimensional tests of both possibilities are 
shown in Figure 3, for three different wave amplitudes. For a wave with maximum amplitude less 
than 3 MPa there is fairly good agreement, however, as the amplitude is increased, as is required 
for ESWT, the linear and quadratic equations of state do not capture the correct behavior. Thus 
we used the full Tait EOS in the fluid domain. 

2.3 Elasticity equations 

In the current work we model bone as a linear isotropic solid. We use the equations (1) together 
with Hooke's law 
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where the spatially- varying scalar coefficients Cij{x,y,z) are determined by the properties of the 
material being modeled. The parameters used for the bone model were found in [42]. 

For an isotropic material we can relate the Cij above to the two Lame parameters, A and /i, that 
are used to model different elastic materials. Ca = A + 2/ifori = 1,...,3, Ca = 2/i for z = 4, . . . , 6, 
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and Cij = A for i 7^ j. Here /i is the shear modulus and A + 2/i is the bulk modulus of the material. 
Note that the A here is different from the A* used to denote the eigenvalues elsewhere in the paper. 

Linear elasticity has been used extensively in the literature to model both trabecular and cortical 
bone [33, 24]. Linear viscoelastic models have also been used for ultrasound wave propagation in 
bone [25]. Our model could be extended to orthotropic models, requiring 9 material parameters, as 
has also been used for bone modeling, e.g. [57]. 



3 Eigenstructure of the hyperbolic system 

The full three-dimensional system of equations (1) models both the nonlinear fluid and the linear 
elastic bone as described in the preceeding sections. This system can be written in quasi- linear form 

qt + X, y, z)q:, + B{q, x, y, z)qy + x, y, z)q^ = 0, (19) 

where A^B and C are the Jacobians of the flux functions in the x, y and z directions respectively. 
For the multi-dimensional methods implemented in Clawpack, we need the solution to the Riemann 
problem along slices in each coordinate direction. Here we provide the details for the solution in the 
x-direction, but the solution in the y and z directions are similar with appropriate perturbations to 
the B and C matrices. 

The corresponding Jacobian for this system in the x-direction is 
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where cr^ss, for example, denotes the partial derivative of a^^ with respect to e^^. In the linear elastic 
case this is simply the coeflacient C13, but the above form also applies to the nonlinear compressible 
equations. The spatial variation in /(g, x, z) and the Jacobian A result from allowing the material 
parameters such as density and elastic moduli to vary in space. The Jacobians in the y- and z- 
directions are similar with the entries perturbed appropriately. 
The eigenvalues for system (20) are 
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When modeling a fluid where the shear stress is zero, there are seven zero- speed eigenvalues since 
o'li2 = cF^is = 0. Only the compressional waves corresponding to A^'^ propagate with nonzero speed. 
Note that with the Tait equation of state (10), 
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In the small amplitude acoustic limit e ^ 0, from (21) we obtain the wave speeds 
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which are the expected waves speeds for compressional waves in the Lagrangian form with this 
equation of state. 

For the elastic sohd, on the other hand, waves 1 and 2 correspond to P- waves while waves 4-6 
correspond to S-waves, and the expected wave speeds are recovered based on the elastic coefficients 
given in Section 2.3. For example, in the x-direction the P-wave speeds are 
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and the S-wave speeds are 
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The corresponding eigenvectors for system (20) are 
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for the P- waves and S-waves, and 
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for the stationary waves. 



3.1 Axisymmetric form of the equations 

We used the two-dimensional axisymmetric form of the equations to generate an initial condition 
for our three-dimensional calculations, as well as for validation of our model. 
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The three-dimensional equations in cyhndrical coordinates are 
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It is interesting to note here that the strain in the 00 direction is non-zero and in this case is cahed 
the hoop strain. A uniform radial displacement is not a rigid body motion, as it would be in the 
two-dimensional plane strain case, but instead produces a circumferential strain. This is because 
the original circumference of the cylinder is 27rr, but when there is a strain in the radial direction 
the circumference grows to 27r(r + i^r), inducing a strain 2™^/27rr = Ur/r. 
The Jacobian for system (29) in the z-direction is 
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(30) 



and has an eigen-structure that is equivalent to the two-dimensional elasticity equations, with the 
addition of a second zero-speed eigenvalue. 
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These equations have the structure 

qt^f{q)r^9{q)z = S{q,r), (31) 

with source terms 

je _ ^ 
r 

put = (32) 

r 

arz 

pwt = . 

r 

In Clawpack, we solve these equations with a fractional-step method. The full problem is split into 
two subproblems that are solved independently. We first solve the homogeneous system obtained by 
setting S' = in (31) using the wave propagation algorithm described in Section 4, and then solve 

qt = S{q,r), (33) 

with an appropriate ODE solver. For (32), we use forward Euler. 



4 Numerical methodology 

We used the wave-propagation algorithms described in [40] and implemented in Clawpack [19] to 
solve the hyperbolic systems of PDEs derived in the preceeding sections. In this section we provide 
the basic details of the numerical methodology and the approximate solution to the Riemann problem 
with a spatially- varying flux function, similar to what was done in [41] We also discuss computational 
issues that require the use of adaptive mesh refinement. 

4.1 Riemann solvers and wave-propagation algorithms 

Recall that the "Riemann problem" is the initial value problem for a 1-dimensional hyperbolic system 
of the form 

gt + /(g,^). = 0, (34) 
with special initial data consisting of two constant states separated by a discontinuity 

/ N f Qi if < 

«0(^^ = | Qr if^>0 • (3^) 

If the flux function is spatially varying then we also use a piecewise-defined flux function with 

The Riemann problem plays a fundamental role in the theory and computation of hyperbolic prob- 
lems, since the Riemann solution consists of waves propagating at constant speeds and can generally 
be computed. For nonlinear systems of equations this is often replaced by an approximate Riemann 
solver as will be discussed below. 

For a linear system of equations qt -\- A{x)qx = the Riemann solution is easily computed in 
terms of the eigenvectors and eigenvalues of the matrices Ai to the left of the interface and Ar to 
the right of the interface. We begin by discussing the linear case with a constant matrix A and turn 
to the variable- coefficient (heterogeneous media) case in Section 4.3. We assume the matrix A is 
diagonalizable, 

A = RAR-^ (37) 
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where R is the matrix of eigenvectors and A is the diagonal matrix of eigenvalues. The Riemann 
solution is computed by decomposing AQ = — Qi as a linear combination of eigenvectors of A, 

m 

AQ = "^aPrP, where a = R-^AQ. (38) 
p=i 

We denote the pth wave by Wp = a^r^^ where p = 1, 2, . . . , m and the number of waves m is 
equal to the number of equations in the system. 

We use finite volume methods in which Qf represents a cell average of the vector q in cell i at 
time tn (still in one space dimension). In Godunov's method the cell average is updated by the 
waves entering the cell from the interfaces to the left and the right, and each wave updates the cell 
average by W^, the jump in q across the wave, multiplied by the distance the wave propagates over 
the time step and divided by the length of the cell, i.e., the cell average is updated by ^^^W^- To 
express the total update to a cell, it is convenient to define matrices and A~ via 

A^ = RA^R-\ where = diag(A^) (39) 

with A+ = max(A, 0) and A~ = min(A, 0). Then the cell average is updated by 

Q7^' = Q7- |^{^+AQi_i/2 + ^-AQi+1/2). (40) 

Here A(5^_i/2 = Qi ~ Qi-i is the jump across the interface at z — 1/2, for example. For a linear 
system this is a generalization of the upwind method and is first order accurate. 
Second order accuracy is achieved by adding in correction fluxes: 

Q^' = Q2- ^(^^^Q + ^~^Q) - - F,_i/2) (41) 

where 



^-1/2 = ^ ( 1 



2 



XPAt 



|A^|Wtv2 (42) 



These terms convert the upwind method into a method of Lax-Wendroff type, matching terms 
through At'^A'^qxx in the Taylor series expansion of the solution at the end of the time step. This 
method generates dispersive errors, however, that can create large nonphysical oscillations near steep 
gradients or discontinuities in a solution, such as shock waves J?o turn this into a "high-resolution" 
method, we use a wave limiter, replacing Wf_i/2 (^^) ^ limited version of the wave. 

The wave Wf_i/2 compared to the corresponding wave from the neighboring Riemann problem, 
either Wf_3^2 if > or Wf_^^^2 if < 0. If the waves are of comparable magnitude the full 
correction term is used for accuracy, but if there is a large discrepancy then the solution is not smooth 
at this point and a limited version is applied. A wide variety of limiters have been developed, and 
we generally use the MC limiter. See [38] or Chapter 6 of [40] for more complete details. 

In two or three space dimensions the idea is the same, but now a 1-dimensional Riemann problem 
must be solved normal to each edge or face of the cell. The resulting waves update the cell averages 
and correction fluxes analogous to (42) are used along with limiters in each direction. 

In addition, to achieve second-order accuracy and good stability properties, it is also necessary 
to use "transverse Riemann solvers" that further modify the correction fluxes F at each cell edge. 
The method described above is based on propagating waves normal to each interface. In reality, the 
waves will propagate in a multidimensional manner and aff"ect cell averages in cells above and below 
those that are directly adjacent to the interface. 

In two space dimensions, each "fluctuation" such as A~ AQi_i^2,j ^Qi-i/2,j that result 

from solving a Riemann problem in the x-direction is split into two pieces using the eigenstructure 
of the coefficient matrix B in the ^/-direction, e.g.. 



(43) 
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These two pieces will modify the correction flux at the edges (i, j — 1/2) and (i, j + 1/2) respec- 
tively to capture the transverse motion of the right-going wave. Similarly, after solving a normal 
Riemann problem in the ^/-direction using the B matrix, transverse problems are solved based on 
the eigenstructure of A. The net effect of all these corrections is to incorporate terms modeling 
the cross- derivative terms BAq^y and ABqyx of the Taylor series expansion in a properly upwinded 
manner. More details can be found in [38] or Chapter 21 of [40]. The transverse correction terms 
are needed for accuracy, but also have the effect of improving the stability limit, allowing a Courant 
number near 1 to be used, relative to the maximum wave speed in any direction. 

In three space dimensions there are two transverse directions for each normal Riemann solve, and 
terms modeling CAq^z^ etc. must also be included. Moreover, "double transverse" terms must be 
included, splitting the result of a transverse solve into eigenvectors of the remaining coefficeient ma- 
trix, and modeling terms such as BCAq^zy The details are presented in [36] and fully implemented 
in Clawpack. 



4.2 The nonlinear fluid Riemann solver 

The compressible fluid equations in Lagrangian form discussed in Section 2.2 can be reduced to 
the quasilinear form (19) in which the Jacobian matrices depend only on q (for a spatially uniform 
fluid) . To apply the wave-propagation algorithm we need to solve the Riemann problem orthogonal 
to each cell interface. For nonlinear problems this is usually done using an approximate Riemann 
solver, e.g., by replacing f{q)x by Aq^^ where the matrix A at each cell interface is chosen based on 
the data Qi and Qr to the left and right. We use the f-wave formulation of the wave-propagation 
algorithm [4], in which the jump in flux f{Qr) — f{Qi) is split into eigenvectors of an approximate 
Jacobian matrix, rather than the jump in Q. This leads to an algorithm that is conservative for any 
choice of approximate Jacobian and also extends naturally to the case of spatially varying fluxes, as 
required near the fluid-solid boundary and discussed further below. 

Rather than choose an approximate Jacobian A and then determining its eigenvectors and eigen- 
values, we simply choose the set of eigenvectors and associated wave speeds based on the data and 
wave forms expected to result from this data. These vectors form a matrix R and we then solve 
up = f{Qr) — f{Qi) for the vector of wave strengths /3. The choice of vectors in R and associated 
wave speeds A implicitly defines the Jacobian approximation A = RAR~^, but this matrix is never 
needed. 

The eigenvectors are taken to be the vectors displayed in (26) and (27). Recall that in the fiuid 
case there are only two nonzero eigenvalues corresponding to the first two eigenvectors. For the 
eigenvector corresponding to < we use A^ = — cr^A evaluated in the left state Q/, while the 
eigenvector corresponding to A^ > is determined using A^ = cr^i\ evaluated in the right state Qr- 
These vectors have nonzero components only in positions 1 and 7 and so the values of and 
can be determined by solving a 2 x 2 system: 



1 1 




in' 1 




r 1 


PiXj pAl 











(44) 



The solution is 



/3^ = 



A/7 -pA? A/1 



PrXl - Pl\] ' pAl - plX] 

The remaining waves do not propagate and do not come into the wave-propagation algorithm. 



(45) 



4.3 The linear elastic Riemann solver 

In the linear elastic material modeling bone, we take a similar approach and again use the f-wave 
formulation of the algorithm. In this case there are six waves with nonzero wave speeds given by the 
eigenvectors in (26). The eigenvectors are independent of q in the linear case, but can be spatially 
varying to represent varying bone structure, so the coefficients Cij in (12) can vary from one grid 
cell to the next. Similar to the nonlinear case described above, to compute the decomposition of 
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Adaptive Grids for ESWT Reflector single Grid vs AIVIR witli 3 levels of Refinement 




200 



Figure 4: a) Resolution of the ellipsoid reflector with different levels of AMR. b) 2d axisymmetric 
calculation. The second is a comparison of the waveform obtained using AMR and a uniform grid. 
The finest grid resolution in the AMR calculation is the same as the resolution on the uniform grid. 



the flux difference into propagating waves we define the three left-going eigenvectors r^'^'^ (with 
the minus sign in (26)) based on the coefficients in the left state, while the right-going eigenvectors 
^2,4,6 defined using the coefficients in the right state. Note that the flux vector f{q) from 
(2), and hence any jump in flux, has zeros in three components which are easily seen to lead to 
= = when the flux difference is written as a linear combination of the eigenvectors, and 
the six remaining components of the flux difference uniquely define the coefficients through 
for the six propagating waves. and are the same as (45), the other weights are 

3 ^ PrKM'-M' .4 ^ A/^-pAfA/^ .... 



4.4 Interfaces and the Cartesian Grid 

In ESWT the pressure wave must propagate through a variety of materials, and in general the 
interfaces between different materials do not align with the grid. In our calculations we use a 
Cartesian grid. To handle grid cells that contain two materials, we perform a weighted average of 
the material properties. The stress-strain relationship in the averaged grid cells is taken to be that 
from linear elasticity, even if one of the materials is fluid. This approach is feasible because we use 
AMR to refine around the interfaces between the two materials. By using a fine enough grid, we are 
able to reduce the error introduced by the weighted average approximation. Figure 4 a) illustrates 
the interface between the fluid and the brass reflector from an axisymmetric calculation. Three grid 
resolutions are shown in this figure: a coarse grid on the right, a level 2 grid that is refined by a 
factor of 4 in each direction in the middle, and the finest grid on the left, where the grid lines are 
not drawn. 

Figure 4 b) shows a comparison between the pressure wave measured at F2 for an AMR calcula- 
tion versus a single grid calculation. The single grid calculation took 269 minutes to complete, just 
over 6 times as long as long as the run using AMR which finished in 44 minutes. These calculations 
were performed serially on a 2.8 GHz dual core AMD Opteron machine with 32 GB of memory. It's 
clear that the two calculations yield comparable pressure waves. The biggest difference is is in the 
direct wave arriving around t = 150, which is not being resolved in the AMR calculation because 
we have refined only in the vicinity of the reflected wave of primary interest. 
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4.5 Adaptive Mesh Refinement 

The pressure waveform found in ESWT contains a very thin region of high pressure that can not be 
resolved without a highly refined mesh. In Figure 4.5 we investigated the effect of grid refinement on 
the shock wave profile and found that with grid resolution greater than 0.25mm, the wave form at 
F2 was not a shock. Note that near the shock we only expect our method to be first order, but the 
solution does converge to a shock as the grid is refined. Our calculations are done with the adaptive 
mesh refinement (AMR) in the style of Collela-Berger-Oliger [6], [5]. The AMR algorithms used in 
Clawpack are more fully described in [7]. For the three-dimensional calculations, a similar AMR 
algorithm is used, as implemented in Chombo. Here we only briefiy review the main ideas. 

The computational domain is covered by a rectangular level 1 grid, typically at a coarse resolution. 
Rectangular patches of the grid may be covered by level 2 grids, refined by some specified refinement 
ratio in each direction. Since we use explicit methods, the CFL condition generally requires that 
the time step be refined by the same factor on the level 2 grids, so several time steps must be taken 
on each level 2 grid for each time step on the level 1 grid. The level 1 grid is advanced first, and 
for each time step on the level 2 grid, ghost cell values around the boundary are filled in either by 
copying from adjacent grids at the same level, or using space-time interpolation from the level 1 grid 
for ghost cells that do not lie in an adjacent grid. This entire procedure is repeated recursively to 
obtain higher levels of refinement, e.g. some portion of the collection of level 2 grids may be covered 
by level 3 grids and so on. 

In order to adaptively refine the grid, it is important to specify appropriate refinement criteria. 
The perturbations to the strain are small, so gradients in the strain are too small to use as reliable 
refinement criteria. However, the small strains result in large changes in the pressure, so we refine 
in the area near the pressure wave. In order to handle the interfaces between two materials, we 
also use large gradients in background density as a secondary refinement criterion. Cells that are 
fiagged as needing refinement are clustered into rectangular patches using the algorithm of Berger 
and Rigoutsis [8] . Regridding is done every few steps on each grid level in order to track propagating 
waves. Regions are automatically de-refined once the wave passes by, since in these regions cells are 
no longer fiagged as needing refinement. 

40 
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Figure 5: Effect of grid size on shock wave profile. As the grid is refined for the same initial 
condition, the shock wave profile steepens. The solution eventually converges to a profile with the 
same magnitude, though the convergence rate is only first order near a discontinuity. 

Figure 4.5 illustrates the behavior of the ESWT waveform as the grid is refined. What is evident 
from these experiments is that a coarse grid will not effectively capture the development of the 
shock, so around the propagating wave, we need at least Ax = .25mm resolution. As the wave 
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steepens into a shock, we no longer expect second order convergence, because in the region around a 
discontinuity, our methods are first order. However, since the discontinuities occur in a smaU region 
of the domain, the overall methodology is still second order. 

In order to efficiently calculate a reasonable ESWT waveform in three dimensions, we utilized 
ChomboClaw [10], which uses the adaptive mesh refinement routines of CHOMBO with the wave 
propagation solvers of Clawpack. This code can be run in parallel using MPI on an NSF TeraGrid 
computer at TACC. We found that the ChomboClaw code scaled up to 128 processors: the same 
experiment ran approximately 4 times faster on 128 processors than on 32. However, we need to 
do further testing on larger systems in order to validate that ChomboClaw will scale as well as 
CHOMBO, which is reported to scale well to 10,000 processors [16]. 

5 Results 

We have used the approach described above to model ESWT pressure waves interacting with three- 
dimensional bone geometries comprised of idealized materials. We have modeled both simple objects 
that have been used in laboratory experiments [34] as well as complex three-dimensional geome- 
tries extracted from CT scans of patient data [21]. Here we present results that demonstrate the 
efficacy of the Lagrangian formulation, as well as examples of calculations performed using real 
three-dimensional geometries. 

The calculations were initialized using pressure data obtained from a 2D axisymmetric calculation 
where we modeled the full geometry of the ellipsoidal reffector. The reffector was modeled using 
linear elasticity with material properties that can be found in [20] . We assumed the ffuid was water 
with the corresponding parameters for the Tait equation of state found in Section 2.1. We saved 
the data at t = 116/is and used this to restart future calculations. For the three-dimensional initial 
condition, we rotated the 2D data about the x-axis. The material properties of averaged bone were 
obtained from [42] and used in the heterotopic ossification, cylinder and sphere. 

We have found in our experiments that interfaces between materials with large impedance dif- 
ferences have the most significant effect on maximum stress and energy deposition. 

5.1 Reflection and focusing 

In Figure 6, we show an axisymmetric calculation of the ESWT wave propagation and focusing in 
water alone, in a domain bounded by the ellipsoidal reflector of the Dornier HM3. Figure 6 a) shows 
the initial spherical propagation of the pressure wave, as well as the grids where the calculation is 
being reflned. The grid must be reflned around the pressure wave as well as the reflector. Figures 6 
b) and c) show the propagation of the wave and evolution of the the adaptive grid structures. At 
later times the grid is only being refined near the pressure wave. The sharp results and absence 
of spurious oscillations in the pressure measurement at F2 indicate that AMR together with our 
Cartesian grid approach enables us to capture the reflection at the interface. 
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Pressure at F2, A x = .25mm 
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Stress at time t = 180 microsec 
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Figure 6: Axisymmetric calculation of the pressure pulse generated by a spherical high-pressure 
bubble centered z = —115 (the focus Fl of the ellipsoidal reflector). Three levels of AMR are 
used and grid lines are shown only on levels 1 and 2. The level 3 grid has a resolution of A2: = 
Ar = 0.25mm. a) At t = 10 the pulse has nearly reached the reflector, b) At t = 70 the incident, 
transmitted, and reflected pulses are visible, c) At t = 180 the reflected pulse has focussed near 
2; = 115 (the focus F2). d) The time history of the pressure at F2. The direct (unreflected) wave 
passes F2 at t ~ 150 and the focussed pulse arrives at t ~ 180. 
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5.2 Axisymmetric sphere 

We used an axisymmetric test problem in order to compare the solutions obtained with the two- 
dimensional and three-dimensional codes. The initial condition for this experiment was an analytic 
form for an ESWT pressure wave used in [54]. In the two-dimensional case, we specified the pressure 
as a function of the radial distance from Fl (-115,0). In the three-dimensional case, we rotated the 
same two-dimensional initial condition about the x-axis. The grid resolution was Ax = .25mm. 

Results with contour lines are in Figure 7. The maximum values in each of the three cases are 
nearly the same, but there are slight discrepancies in the contour lines. Figure 8 shows a ID slice 
of the maximum shear calculation in the 2D and 3D codes, which makes it clear that the peak of 
maximal shear stress is in the same location and has the same value. The general shape of the 
maximum stress deposition pattern are similar in both cases. The difference in the two solutions 
is likely caused by the solid wall boundary condition that is used at r = 0. Only waves that are 
propagating normal to that boundary are perfectly refiected, otherwise some error is generated. 
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Figure 7: Results from calculation of a shockwave interacting with an acrylic sphere. The left column 
shows 2D axisymmetric results and the right column shows a corresponding cross section of full 3D 
calculation. Top: Maximum compression. Middle: Maximum tension. Bottom: Maximum shear. 
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Slice of Maximum Sliear Stress in Spliere along x,y=0 
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Figure 8: Comparison of maximum shear stress from 2D and 3D calculations as a function of x 
along y = z = {). The difference in the results is likely caused by averaging of the initial condition 
onto the 3D domain and the boundary conditions on the axisymmetric calculation at r = 0, but the 
two calculations predict comparable location and magnitude of maximal shear stress deposition. 



5.3 Nonunions 

ESWT has recently been used in the treatment of non-unions or bone fractures that fail to heal [9]. 
One question that is of interest to clinicians is whether or not the angle of treatment has an effect on 
healing. We assume that healing is related to the magnitude of stress applied near the treatment area, 
although the connection between the applied force and biological response is not yet understood. In 
the fluid there is no shear stress. However, at the liquid-solid interfaces, shear stresses are generated 
by the shockwave and stimulate motion both at the surface and within the material. The motion 
of the biological materials (e.g. the periosteum, interstitial fluid, mechanotransduction) is likely to 
be important in the healing process [28, 49, 62, 13, 30, 48, 58], and modeling the magnitude and 
location of the stress deposition is a good first step toward understanding the shear and tensile 
displacements caused by ESWT. We should stress, however, that the healing mechanisms are not 
well understood and we are not claiming that magnitude of the applied stress is the most important 
or only biological mechanism involved in the healing process. As mentioned in Section 1, several 
studies have indicated that cyclic application of mechanical loading leads to the generation of new 
bone. The work of Isaksson, et. al. [30], indicates that the most accurate predictors for bone healing 
are those based on shear strain and fluid flow, however, there is no single model that can predict all 
features of the healing process, so more work is necessary [44]. 

In an actual treatment, the clinician generally sets up the device so that the focus is aligned with 
the ailment. For example, in the case of a broken bone, the clinician will set up the device so that 
F2 is in the center of the break. However, given the heterogeneous media, it is not clear that the 
maximal stresses will actually be observed at F2, as would be expected in pure water. We used our 
model to investigate the location of maximal stress deposition relative to F2. In these calculations 
we considered two different geometries, a complete cylinder, representing the long shaft of a healthy 
bone, and a broken cylinder, representing a nonunion. The results from calculations where the 
idealized bone was perpendicular to the direction of the pressure wave front are shown in figures 9 
and 10. We found that the break has a significant impact on the location of stress deposition. 
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We used these geometries to perform a variety of experiments. We rotated the direction of 
treatment by 45 and 60 degrees relative to the x-axis and calculated both the magnitude of the 
maximum compressive, tensile and shear stresses, as well as the distance from the focus F2 of the 
device. 

In the case of the broken cylinder, the maximum stress deposition in the direct experiment is 
similar to that of the unbroken cylinder, except that the there are two locations of maximal stress 
deposition on either side of the break. The pressures in the bone are larger than in the fluid due to 
reflection at the fluid-solid interface, so the contours of maximum stress are concentrated on either 
side of the gap. The location along the x-axis is nearly the same as in the unbroken cylinder, and 
the distances from the ideal focal point, F2, are also similar. 

As the angle of treatment is varied, there is less of a shift in the z-direction for the shear 
and compressive stresses. This is caused by the impedance diff'erence between the fluid and solid 
material at the gap, which is located close to F2. If the gap were shifted along the z-axis from the 
focal point, there would be a corresponding shift in the location of maximum shear and compression. 
Geometrically, the shape of the regions of compressive and shear stress are quite diff'erent from the 
direct case. Instead of being an ellipsoidal shape, the regions are compressed into the corner of the 
lower-half of the cylinder. Again, this is caused by the impedance jump at the fluid-solid interface. 
The region of maximum tension deposition is similar to that of the unbroken cylinder case, though 
it is also affected by the gap and the tension is concentrated on the upper half of the cylinder. 

It is clear from the literature [48, 49, 18, 13, 30, 11, 26], that mechanical loading is important in 
bone healing. The implication of our computational experiments is that the angle of treatment will 
affect stress deposition and therefore may be important in treatment optimization. For example, in 
order to maximize shear stress at the tissue-bone interface, our preliminary computations indicate 
that it might be best to treat the patient at an oblique angle. However, if the goal is to maximize 
shear stress in the gap of the broken bone, then treating the patient at a 90 degree angle may be 
better than treating at either the 45 or 60 degree angle. We stress however that the biological 
mechanisms must be better understood and more experiments must be done in conjunction with 
laboratory and clinical treatments before these calculations could be used to make specific clinical 
recommendations . 

In Figure 11 we show two-dimensional slices of a calculation with a more realistic, but still 
idealized, long bone geometry. In a the shaft of a long bone, there is a marrow-filled canal running 
through the center. Marrow is typically modeled as a viscoelastic material [42], but for this first 
approximation we just used a fluid-filled canal. The impedance difference in the two materials is 
similar and therefore illustrates the behavior that we are interested in, the change in maximal stress 
deposition. In contrast to the solid cylinder above, the contours of maximal stress are concentrated 
in the front side of the hollow cylinder. Figure 11 b) shows that there are also two regions of 
additional stress concentration in the backside of the hollow cylinder. This example highlights the 
importance in understanding where these impedance jumps occur in order to optimally treat the 
patient. 

5.4 Heterotopic Ossification 

A heterotopic ossiflcation (HO) is a growth of bone-like material in soft tissue. HOs often grow 
spontaneously in tissue that has been traumatized due to injury or amputation. An example of an 
HO is shown in Figure 12 a), which shows the pelvis and an HO, using data extracted from a patient's 
CT scan. In this case, the HO has grown around the right hip joint and is inhibiting the patient's 
range of motion. The goal of the HO treatment is not to pulverize the ossiflcation, but to break 
up the adhesion between the HO and the joint, in order to restore the patient's range of motion. 
There is no clear division between the HO and bone in the CT scan because both are composed 
of materials that have similar densities. However, the HO does not have the same woven structure 
that is present in bone, so the two will likely have different material properties, even though the 
densities are similar. This similarity means that we are uncertain as to how strong the connection or 
adhesion is between the HO and the bone, which will directly impact the number of shocks needed 
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Figure 9: Three-dimensional results for the direct treatment of a complete cylinder. This figure 
shows 2D slices of maximum compression, tension and shear along y=0 for treatment where the 
ESWT wave propagates along the x-axis, as indicated by the arrow. The dot illustrates the location 
of F2. 
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Figure 10: Three-dimensional results for the direct treatment of a broken cylinder. This figure shows 
2D slices of maximum compression along y=0 for treatment along the x-axis, 45 degree rotation and 
60 degree rotation about the y-axis. The arrows indicate the angle of treatment in each case and 
the dot illustrates the location of F2. 



to restore the patient's range of motion. 

We are able to use our model to investigate the effect of the angle of treatment on the observed 
stresses in the region near the HO. Since the composition and material properties of the ossification 
are not well understood, we can also use the model to vary the material properties of the ossification 
and investigate the sensitivity of the results to these parameters. We found that both the strength 
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Figure 11: An additional interface calculation showing the more realistic treatment of a cylinder 
with a fluid-filled cavity, a) is a slice along z=0 showing the concentration of stress in the front of 
the idealized bone, with additional smaller pockets of maximum stress due to refiection in the back 
half of the bone, b) is a slice along y=0 which further demonstrates the stress concentration in the 
first half of the bone. The arrow in b indicates the direction of ESWT wave propagation and the 
dot indicates the location of F2. 



of the connection between the HO and bone, as well as the composition of the HO, had a significant 
effect on the location of maximum stress in the object [21]. 

It is challenging to infer anything meaningful from the images in the full three-dimensional 
calculation in Figure 12 a), so we have also included a two-dimensional slice of the maximum shear 
in Figure 12 b). Here the gray regions represent the bone-like HO material and we assume any gaps 
are filled with fiuid. It is clear that the interior of the ossification is complex and contains may 
fiuid-filled pockets that affect, in this case, the location of the maximum shear stress. 

Given the complex nature of the HO and subsequent difficulty interpreting the three-dimensional 
results, we have also used an idealized ossification to investigate some facets of the treatment. One 
example of this is shown in Figure 13, where we have simulated a case where the ossification (the 
crescent in the two-dimensional images) is not strongly attached to the bone (the cylinder) and 
calculated the maximal shear stress as a result of two different treatment angles. Figure 13 a) 
illustrates the result when the ESWT device is aimed orthogonal to the gap between the HO and 
the cylinder. Figure 13 b) is the result when the device is aimed so the Shockwaves propagate parallel 
to the gap between the HO and bone. It has been indicated that maximum shear stress is important 
in causing the HO to break, so it is desirable to deposit the maximum amount of shear as close to 
the HO/bone interface as possible. In this case, it is better to treat the HO in the direction indicated 
in Figure 13 b), since the shear stress is concentrated along the gap. 

According to our computational results, the pockets of fluid within an HO and strength of 
adhesion to the bone surface will affect the stress deposition and therefore the location of the 
eventual break in the ossiflcation. Further investigation is required to be conclusive, but our results 
indicate that if the fluid pockets are in the propagation path of the shock wave, they may cause 
the maximum stresses to occur away from the adhesion site, making the treatment less effective. 
In reality, the composition of the HO is unknown and we do not have a good characterization for 
the material properties of the ossiflcations. However, the strong impedance mismatch between fluid 
and bone, as well as the inability of the fluid to support shear stress, indicate that the presence 
of fluid-fllled pockets will have an effect on the stress deposition. We should note here that our 
modeling work does not take into account the propagation of successive shocks or failure models 
within the material, which should ultimately be incorporated in order to the determine the optimal 
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treatment. This is an area for future work. 




Figure 12: a) The 3D CT patient data ihustrating the heterotopic ossification (blue) attached to 
the right hip joint (green), b) A shce at x=115 of the 2D calculation shows how the pockets of fluid 
lead to stress concentration in the substructure of the ossification, the dot indicates the location of 
F2 and the direction of treatment is into the page. 




Figure 13: Calculations for an idealized ossification that demonstrates the difference in shear stress 
deposition when treating the HO from different directions. Since the goal is to disrupt the HO at 
the interface between the bone and the HO, b) indicates that it would be optimal send shock waves 
parallel to the break, instead of perpendicular to it. The arrows indicate the direction of ESWT 
propagation and the dots indicate the location of F2. 



6 Conclusion 

In this paper we have proposed a new model for ESWT. We have demonstrated that the Tait 
equation of state is sufficient for the pressures that arise in ESWT. We have shown that the fluid 
and solid can be modeled with the same set of Lagrangian equations since the particle displacements 
are small. This approach allowed us to utilized existing numerical methodology, consisting of high- 
resolution shock-capturing methods together with Adaptive Mesh Refinement, to efficiently calculate 
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solutions to these equations for a variety of idealized biological problems. We have also demonstrated 
that we can effectively handle interfaces between different materials on Cartesian grids. Using 
this methodology we were able to explore, even in geometrically complicated structures, how the 
interfaces between the fluid and solid materials affect the distribution of maximal stress in several 
problems of clinical interest. 

Maximizing stress in specific regions seems important in both the healing and destruction of 
biological tissues. Shear stress is thought to be play a role in the stimulation of biological tissues 
[60, 28, 23, 26, 13, 18, 48]. Mechanical loading is thought to play a role in the formation of bone tissue, 
and as discussed in Section 1, shear and compressive displacements generated by loading influence 
bone healing [48, 49, 18, 13, 30, 11, 26]. Shear stress is also important in predicting the break up 
of kidney stones [54] . Tension leads to the formation of cavitation bubbles that lead to damage and 
crack formation in kidney stones when subjected to shock wave lithotripsy. These cracks are thought 
to be one mechanism that initiates kidney stone comminution, or the pulverization of the stone. The 
role of cavitation in the treatment of heterotopic ossifications is not yet well understood, but given 
the similarity of HO treatments to lithotripsy, it is reasonable to hypothesize that cavitation will 
impact the treatment. 

The model we have developed has been used to investigate idealized non-unions and heterotopic 
ossifications, and we have shown a few examples to illustrate this. A broader range of calculations 
are available in [20]. These specific examples are being studied in more detail and papers with 
additional work on these problems are now in preparation [21, 34]. 

The focus of this paper has been the effect that material interfaces between tissue and bone have 
on the transmission, reflection, and focusing of the shock wave. Very simple models have been used 
for the material on each side of the interface: compressible fluid with a Tait equation of state in the 
tissue and linear isotropic elasticity in the bone. We believe that this level of macroscopic modeling 
can already reveal interesting features of the stress that may be clinically important. In particular, 
focusing may occur in regions displaced from where it would be observed in pure water, and mode 
conversion at an interface can generate shear waves in the bone that are not present in the focusing 
shock wave in fluid. 

To consider the effect of stress on individual osteocytes, a much more detailed model would 
be necessary that is beyond the scope of this work. In particular, this would require modeling the 
microscale fluid-filled canaliculi within the bone through which the osteocyte processes extend. Work 
is currently underway in this direction, and also on intermediate levels of modeling in which the bone 
is modeled as an orthotropic poroelastic material. These equations can be solved with essentially 
the same high resolution finite volume methods used here, after implementing a more complicated 
Riemann solver [37], and with the same software for adaptive mesh refinement. Another possible 
extension is to investigate viscoelastic tissue models that may be superior to the Tait equation for 
water that is currently used. 
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